The anisotropy of granular materials 
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The effect of the anisotropy on the elastoplastic response of two dimensional packed samples of 
polygons is investigated here, using molecular dynamics simulation. We show a correlation between 
fabric coefficients, characterizing the anisotropy of the granular skeleton, and the anisotropy of the 
elastic response. We also study the anisotropy induced by shearing on the subnetwork of the sliding 
contacts. This anisotropy provides an explanation to some features of the plastic deformation of 
granular media. 
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I. INTRODUCTION 

The mechanical behavior of granular materials has 
been largely investigated using constitutive models. 
These are empirical relations which are based on, for ex- 
ample, laboratory tests of soil specimens. On the other 
hand, the investigation of the soils at the grain scale, us- 
ing discrete element modeling, has become possible in re- 
cent years. These models have provided a valuable under- 
standing of many micro-mechanical aspects of soil defor- 
mation. In particular, numerical simulations of packings 
of disks evidence that the stress applied on the bound- 
ary of the assembly is transmitted through a heteroge- 
neous network of interparticle contacts [1]. The geomet- 
ric change of this network during deformation shows a 
structural anisotropy induced by shearing [2-4] . Numer- 
ical simulations have also shown a relevant number of 
contacts reaching the sliding condition, even when the 
sample is isotropically compressed [5,1]. However, little 
work has been done in order to connect the sliding at the 
contacts to the plastic deformation of granular materi- 
als. The aim of this paper is to combine the continuous 
and the discrete approaches in the investigation of the 
anisotropic response of granular materials, taking into 
account the anisotropy induced in both sliding and non- 
sliding contacts. 

We will numerically study the elasto-plastic response 
of a two-dimensional granular model material. The in- 
terparticle forces include elasticity, viscous damping and 
friction with the possibility of slippage. The system is 
driven by applying stress controlled tests at the bound- 
ary particles. The results show that the traditional fabric 
tensor is not sufficient to describe the complex elasto- 
plastic response of granular materials. Additional pa- 
rameters taking into account the anisotropy of the sub- 
network of the sliding contacts are necessary to include 
in the description of the overall plastic deformations. 

This paper is organized as follows: The details of the 
particle model are presented in Sec. II. The contact 
forces are implemented by a Coulomb friction criterion, 
and the stress is controlled by the application of suit- 



able forces at the boundary particles. The calculation of 
the constitutive relations is presented in Sec. III. Here 
we discuss the results in the framework of the Drucker- 
Prager theory of elasto-plasticity. In Sec. IV we study 
the effect of the anisotropy of the contact network in the 
incremental elastic response of the assembly. In Sec. V 
some internal variables are introduced in the constitutive 
relations. These variables take into account the effect of 
the anisotropy induced in the subnetwork of the sliding 
contacts on the plastic deformation of the assembly. 



II. MODEL 

We present here an extension of those two-dimensional 
discrete element methods which have been used to model 
granular materials via polygonal particles [6,7]. The de- 
tails of the particle generation, the contact forces, the 
boundary conditions and the molecular dynamics simu- 
lations are presented in this section. 



A. Generation of polygons 

The polygons representing the particles in this model 
are generated by using the method of Voronoi tessella- 
tion [7]. This methods is schematically shown in Fig. 1: 
First, a regular square lattice of side £ is created. Then, 
we choose a random point in each cell of the rectangular 
grid. Then, each polygon is constructed assigning to each 
point that part of the plane that is nearer to it than to 
any other point. The details of the construction of the 
Voronoi cells can be found in the literature [8,9]. 

Using the Euler theorem, it has been shown analyt- 
ically that the mean number of edges of this Voronoi 
construction must be six [9] . The number of edges of the 
polygons is distributed between 4 and 8 for 98.7% of the 
polygons. We also found that the orientational distribu- 
tion of edges is isotropic, and the distribution of areas 
of polygons is symmetric around its mean value £ 2 . The 
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probabilistic distribution of areas follows approximately 
a Gaussian distribution with variance of 0.36£ 2 . 




FIG. 1. Voronoi construction used to generate the convex 
polygons. The dots indicate the point used in the tessellation. 
Periodic boundary conditions were used. 



B. Contact forces 



In order to give a three-dimensional picture of this 
model, one can consider the polygons as a collection of 
prismatic bodies with randomly-shaped polygonal basis. 
We assume that all the bodies have the same thickness 
L. The force between two polygons is written as F — Lf 
and the mass of the polygons is M = Lm. 

In reality, when two elastic bodies come into contact, 
they have a slight deformation in the contact region. In 
the calculation of the contact force we suppose that the 
polygons are rigid, but we allow them to overlap. Then, 
we calculate the force from this virtual overlap. 

The first step for the calculation of the contact force is 
the definition of the line representing the flattened con- 
tact surface between the two bodies in contact. This is 
defined from the contact points resulting from the in- 
tersection of the edges of the overlapping polygons. In 
most cases, we have two contact points, as shown in the 
left of Fig. 2. In such a case, the contact line is de- 
fined by the vector C = C\Ci connecting these two in- 
tersection points. In some pathological cases, the inter- 
section of the polygons leads to four or six contact points. 
In these cases, we define the contact line by the vector 
C = C^cf 2 + C^Ul orC = CYct + C^Ul + C^ct, respec- 
tively. This choice guarantees a continuous change of the 
contact line, and therefore of the contact forces, during 
the evolution of the contact. 
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FIG. 2. Contact points d before (left) and after the for- 
mation of a pathological contact (right). The vector denotes 
the contact line, t represents the time step. 

The contact force is separated as f c = f e + f v , 
where f e and f v are the elastic and viscous contribu- 
tion. The elastic part of the contact force is decomposed 
as f e — fn^ c + ftt c - The calculation of these compo- 
nents is explained below. The unit tangential vector is 
defined as t c — C/\C\, and the normal unit vector n c is 
taken perpendicular to C. The point of application of 
the contact force is taken as the center of mass of the 
overlapping polygons. 

As opposed to the Hertz theory for round contacts, 
there is no exact way to calculate the normal force 
between interacting polygons. An alternative method 
has been proposed in order to model this force [6]. In 
this method, normal elastic force is calculated as = 
—k n A/L c where k n is the normal stiffness, A is the over- 
lapping area and L c is a characteristic length of the poly- 
gon pair. Our choice is L c = \C\. This normalization is 
necessary to reflect the fact that the normal clastic force 
must be proportional to an overlapping length as it was 
shown in bidimensional Hertzian contacts [10]. 

In order to model the quasi-static friction force, we cal- 
culate the elastic tangential force using an extension of 
the method proposed by Cundall-Strack [11]. An elastic 
force / t e = —k t Ax t proportional to the elastic displace- 
ment is included at each contact. k t is the tangential 
stiffness. The elastic displacement Axt is calculated as 
the time integral of the tangential velocity of the contact 
during the time where the elastic condition |/ t e | < /i/^ is 
satisfied. The sliding condition is imposed, keeping this 
force constant when |/ t e | — fif^. The straightforward cal- 
culation of this elastic displacement is given by the time 
integral starting at the beginning of the contact: 



V C t (t')Q([lfe-\f?\)dt', 



(1) 



where O is the Heaviside step function and v% denotes 
the tangential component of the relative velocity v c at 
the contact. 



LUi X li + UJj X £j . 



(2) 
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Here Vi is the linear velocity and Ui is the angular ve- 
locity of the particles in contact. The branch vector ti 
connects the center of mass of particle i with the point of 
application of the contact force. Eq. (1) allows to keep 
Ax t at a length such that |/ t e | agrees with /j,f% during 
the sliding condition. 

Damping forces are included in order to allow for rapid 
relaxation during the preparation of the sample, and to 
reduce the acoustic waves produced during the loading. 
These forces are calculated as 

/^ = -m(7„<n c + m c i c ), (3) 

with m = (1 /mi + l/m^) -1 , the effective mass of the 
polygons in contact. n c and t c are the normal and tan- 
gential unit vectors defined before, and j n and 74 are the 
coefficients of viscosity. These forces introduce time de- 
pendent effects during the cyclic loading. However, we 
will show that these effects can be arbitrarily reduced 
by increasing the loading time, as it corresponds to the 
quasi-static approximation. 

C. Molecular dynamics simulation 

The evolution of the position Xi and the orientation ipi 
of the polygon i is governed by the equations of motion: 

c b 

im = ^xf9 + ^x^. (4) 

c b 

Here m,i and 7j are the mass and moment of inertia of 
the polygon i. The first summation goes over all particles 
in contact with this polygon. According to the previous 
section, these forces f c are given by 

f c = -(k n A/L c + im v c n )n c - (Ax c t + ln mv c t )t c , 

(5) 

The second summation on the left hand side of Eq. (4) 
goes over all the edges of the polygons in contact with 
the external contour of the assembly. We apply a force 

= -a 1 Ax b 2 x\ + a 3 Ax\x2 - m^vl. (6) 

on each selected segment T b — Ax\x\ + Ax b 2 X2 of the 
external contour of the assembly. Here x\ and xi are the 
unit vectors of the Cartesian coordinate system. o\ and 
03 are the components of the stress we want to apply on 
the sample, as we see in Subsec. III. m; and are the 
mass and the velocity of the particle i belonging to the 
boundary, ft is the vector connecting the center of mass 
of the boundary particle i to the center of the boundary 
segment T b . 



We use a fifth-order Gear predictor-corrector method 
for solving the equation of motion [12]. This algorithm 
consists of three steps. The first step predicts position 
and velocity of the particles by means of a Taylor expan- 
sion. The second step calculates the forces as a function 
of the predicted positions and velocities. The third step 
corrects the positions and velocities in order to optimize 
the stability of the algorithm. This method is much more 
efficient than the simple Euler approach or the Rungc- 
Kutta method, especially for problems where very high 
accuracy is a requirement. 

The parameters of the molecular dynamics simulations 
were adjusted according to the following criteria: 1) guar- 
antee the stability of the numerical solution, 2) optimize 
the time of the calculation, and 3) provide a reasonable 
agreement with the experimental data. 

There are many parameters in the molecular dynam- 
ics algorithm. Before choosing them, it is convenient to 
make a dimensional analysis. In this way, we can keep 
the scale invariance of the model and reduce the param- 
eters to a minimum of dimensionless constants. First, 
we introduce the following characteristic times of the 
simulations: the loading time to, the relaxation times 
t n = l/7„, tt — I/74, h — I/76 and the characteristic 
period of oscillation t s = \Jk n / 'pi 2 of the normal con- 
tact. 

Using the Buckingham Pi theorem [13], one can show 
that the strain response, or any other dimensionless vari- 
able measuring the response of the assembly during load- 
ing, depends only on the following dimensionless param- 
eters: ai = t n /t s , a 2 = t t /t s , a 3 = t b /t s , a 4 = to/t s , 
the ratio £ = kt/k n between the stiffnesses, the friction 
coefficient \x and the ratio <Ji/k n between the stresses and 
the normal stiffness. 

The variables Qj will be called control parameters. 
They are chosen in order to satisfy the quasi-static ap- 
proximation, i.e. the response of the system does not 
depend on the loading time, but a change of these param- 
eters within this limit does not affect the strain response. 
ct2 = 0.1 and «2 = 0.5 were taken large enough to have 
a high dissipation, but not too large to keep the numer- 
ical stability of the method. a 3 = 0.001 is chosen by 
optimizing the time of consolidation of the sample dur- 
ing the application of the confining pressure. The ratio 
a a = to/t s = 10000 was chosen large enough in order to 
avoid rate-dependence in the strain response, correspond- 
ing to the quasi-static approximation. Technically, this 
is performed by looking for the value of 04 such that a 
reduction of it by half makes a change of the stress-strain 
relation less than 5%. 

The parameters £ and /i can be considered as material 
parameters. They determine the constitutive response of 
the system, so they must be adjusted with the experi- 
mental data. In this study, we have adjusted them by 
comparing the simulation of biaxial tests of dense polyg- 
onal packings with the corresponding tests with dense 
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Hostun sand [14]. First, the initial Young modulus of 
the material is linearly related to the value of normal 
stiffness of the contact. Thus, k n — WOMPa is chosen 
by fitting the initial slope of the stress-strain relation in 
the biaxial test. Then, the Poisson ratio depends on the 
ratio ( = k t /k n . Our choice k t — 52.8MPa gives an 
initial Poisson ratio of 0.07. This is obtained from the 
initial slope of the curve of volumetric strain versus ax- 
ial strain. The angles of friction and the dilatancy are 
increasing functions of the friction coefficient fx. Taking 
/i = 0.25 yields angles of friction of 30 — 40 degrees and 
dilatancy angles of 20 — 30 degrees. The experimental 
data yields angles of friction between 40 — 45 degrees and 
dilatancy angles between 7—14 degrees. A better adjust- 
ment would be made by including different void ratios in 
the simulations, but this is out of the scope of this work. 



III. STRESS-STRAIN RELATION 

The characterization of the macroscopic state of a 
granular material in static equilibrium is usually given 
by the Cauchy stress tensor. The average of this tensor 
over the assembly leads to [15] 



cm = 



b 



(7) 



The sum goes over all the forces acting over the bound- 
ary of the assembly. x b is the point of application of the 
boundary force f°. This force is defined in Eq. (6). A is 
the area enclosed by the boundary. The sum goes over 
all the boundary forces of the sample. Inserting Eq. (6) 
in Eq. (7) and taking the equilibrium condition Vi = 
leads to 

-<7 1 E& x b^Ub (J 3 J2b x b^ x b /g-, 

E& yb&yb 0-3 E& Vb^xb \ ' 

Those sums can be converted into integrals over closed 
loops and the calculation of such integrals leads to 



1 

A 



a = 



o x 
a 3 



(9) 



Thus, the stress controlled test is restricted to stress 
states without off-diagonal components. So we can sim- 
plify the notation introducing the pressure p and the de- 
viatoric stress q in the components of the stress vector 



(10) 



In the same way, the incremental strain tensor can be 
calculated from the average of the displacement gradient 
over the area of the RVE. It has been shown [16] that it 
can be transformed in a sum over the boundary of the 
sample 
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Here Au b is the displacement of the boundary segment, 
that is calculated from the linear displacement Ax and 
the angular rotation Acj) of the polygons belonging to it, 
according to 

Au b = Ax + A<f x £. 

From the eigenvalues dei de% of de^ 
umetric and deviatoric components of the strain as the 
components of the incremental strain vector: 



(12) 

we define the vol- 



di 



de 




dei + de3 


dj 




de\ — de% 



(13) 



By convention de > corresponds to a compression of 
the sample. We are going to assume a rate-independent 
constitutive relation between the incremental stress and 
incremental strain tensor. In this case the incremental 
relation can generally be written as [17]: 



de = M(§, a)da, 



(14) 



where 8 is the unitary vector defining a specific direction 
in the stress space: 



da 



cos 9 
sin6> 



\da\ = ^dp 2 + dq 2 . (15) 



The constitutive relation results from the calculation 
of de(9), where each value of 8 is related to a particular 
mode of loading. Some special modes are listed in Table 
I. 

In order to compare the resulting incremental response 
to the elasto-plastic theory, it is necessary to assume that 
the incremental strain can be separated into an elastic 
(recoverable) and a plastic (unrecoverable) component: 



de = d~e e + deP, 
dl e = D^i^da, 
de p = J(9, a)dd. 



(16) 
(17) 
(18) 



0° 


isotropic compression 


dp > 


dq = 


45° 


axial loading 


dai > 


da A = 


90° 


pure shear 


dp = 


dq > 


135° 


lateral loading 


dai = 


da 3 > 


180° 


isotropic expansion 


dp < 


dq = 


225° 


axial stretching 


dai < 


da 3 = 


270° 


pure shear 


dp = 


dq < 


315° 


lateral stretching 


dai = 


da 3 < 



TABLE I. Principal modes of loading according to the ori- 
entation of 6 
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Here, D 1 is the inverse of the stiffness tensor D, and 
J = M - D^ 1 the flow rule of plasticity [18]. They can 
be obtained from the calculation of de e (9) and di p (0) as 
we will see below. 

The method presented here to calculate the strain re- 
sponse has been used on sand experiments [19]. It was 
introduced by Bardet [20] in the calculation of the incre- 
mental response using discrete element methods. This 
method will be used to determine the elastic de e and 
plastic de p components of the strain as function of the 
stress state a and the stress direction 9. 

First, it is isotropically compressed until it reaches the 
stress value o\ = 03 = p — q. Next, it is subjected to 
axial loading in order to increase the axial stress o\ to 
p + q. When the stress state with pressure p and devia- 
toric stress q is reached, the sample is allowed to relax. 
Loading the sample from a to a + da the strain incre- 
ment de is obtained. Then the sample is unloaded to a 
and one finds a remaining strain de p , that corresponds 
to the plastic component of the incremental strain. For 
small stress increments the unloaded stress-strain path is 
almost elastic. Thus, the difference de e = de — de p can be 
taken as the elastic component of the strain. This pro- 
cedure is implemented on different " clones" of the same 
sample, choosing different stress directions and the same 
stress amplitude in each one of them. 

This method is based on the assumption that the strain 
response after a reversal loading is completely elastic. 
However, numerical simulations have shown that this as- 
sumption is not strictly true, because sliding contacts are 
always observed during the unload path [21,5]. In order 
to overcome this difficulty, Calvetti et al. [21] calculate 
the elastic part by removing the frictional condition set- 
ting fj, = co, and measuring the purely elastic response 
de ns of the assembly. Then the plastic component of the 
strain can be calculated as de? = de — de ns . 

In our simulations, we have observed that the plastic 
deformation during the reversal stress path is less than 
1% of the corresponding value of the elastic response. 
Within this margin of error, the method proposed by 
Bardet can be taken as a reasonable approximation to 
describe the elasto-plastic response. 

Fig. 3 shows the load-unload stress paths and the cor- 
responding strain response when an initial stress state 
with <j\ — 200k Pa and 03 = 120k Pa is chosen. The end 
of the load paths in the stress space maps into a strain 
envelope response de(9) in the strain space. Likewise, the 
end of the unload paths maps into a plastic envelope re- 
sponse de p (9). The yield direction (f> can be found from 
this response, as the direction in the stress space where 
the plastic response is maximal. In this example, this is 
around 9 = 87.2°. The flow direction ip is given by the 
direction of the maximal plastic response in the strain 
space, which is around 76.7°. The fact that these direc- 
tions do not agree reflects a non- associated flow rule, as 
it is observed in experiments on realistic soils [19]. 
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FIG. 3. Stress - strain relation resulting from the load - un- 
load test. Dotted lines represent the paths in the stress and 
strain spaces. The dash-dot line shows the strain envelope 
response and the solid line is the plastic envelope response. 



IV. ELASTIC RESPONSE 



Hooke's law of elasticity states that the stiffness ten- 
sor of isotropic materials can be written in terms of two 
material parameters, i.e. the Young modulus E and the 
Poisson ratio v [22] However, the isotropy is not fulfilled 
when the sample is subjected to deviatoric loading. In- 
deed, numerical simulations [2,23] and photo-elastic ex- 
periments [24] on granular materials show that loading 
induces a significant deviation from isotropy in the con- 
tact network. 
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A. Anisotropy of the contact network 

The anisotropy of the granular sample can be char- 
acterized by the distribution of the orientations of the 
branch vectors £, as shown in Fig. 4. Each branch vec- 
tor connects the center of mass of the polygon to the 
center of application of the contact force. Fig. 4 shows 
the branch vectors for two different stages of loading. 
The structural changes of micro-contacts are principally 
due to the opening of contacts whose branch vectors are 
oriented nearly perpendicular to the loading direction. 
Let us call Q((p)Aip the number of contacts per parti- 
cle whose branch vector is oriented between the angles ip 
and ip+Aip, measured with respect to the direction along 
which the sample is loaded. The anisotropy of the con- 
tact network can be accurately described by a truncated 
series expansion. 



No 
2tt 



a a + a\ cos(2iy9) + d2 cos(4(^)] . 



(19) 




n(<p) 



FIG. 4. The lines show the branch vectors for 
o-i=o-3 = IGOkPa (left) and a x = 272kPa and 03 = 48kPa 
(right). The orientational distribution of branch vectors is 
shown for both cases. 
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FIG. 5. Fabric coefficients of the distribution of the contact 
normal vectors. They are defined in Eq. (19). 

Here N = Ngao is the average coordination number 
of the polygons, whose initial value Nq = 6.0 reduces as 
the load is increased. The parameters ao, a\ and an can 



be related respectively to the zeroth, second and fourth 
order fabric tensor defined in other works to characterize 
the orientational distribution of the contacts [2,10]. Here, 
they will be called fabric coefficients. The dependence of 
the fabric coefficients on the stress ratio q/p is shown 
in Fig. 5. We observe that for stress states satisfying 
q < OAp there are almost no open contacts. Above this 
limit a significant number of contacts are open, leading 
to an anisotropy in the contact network. Fourth order 
terms in the Fourier expansion are necessary in order to 
accurately describe this distribution. 



B. Anisotropic stiffness 

We will investigate the effect of the anisotropy of the 
contact network on the stiffness of the material. The 
most general linear relation between the incremental 
stress and the incremental elastic strain for anisotropic 
materials is given by 



d<Tij — Dijkide^, 



(20) 



where D^ki is the stiffness tensor [22]. Since the stress 
and the strain are symmetric tensors, one can reduce 
their number of components from 9 to 6, and the num- 
ber of components of the stiffness tensor from 81 to 
36. Further, by transposing Eq. (20) one obtains that 
Dijki — Djuk, which reduces the constants from 36 to 
21. In the particular case of isotropic materials, it has 
been shown that the number of constants can be reduced 
to 2 [22]: 



[(1 - v)daij - v6ijda kk }. 



(21) 



Here E is the Young modulus and v the Poisson ratio. 
The stress-strain relation of Eq. (20) has been inverted 
to compare to the elasto-plastic relation of Eq. (17). The 
description of the general case of the anisotropic elasticity 
with 21 constants does not seem trivial. However, since 
we consider here only plane strain deformations, we can 
perform further simplifications. We take a coordinate 
system oriented in the principal stress-strain directions. 
Thus, the only non-zero components are dan = do~\ and 
do-33 = d<73 for the stress and den = do\ and ^£33 = das 
for the strain. The anisotropic elastic tensor connecting 
these components contains only three independent pa- 
rameters. We can write Eq. (20) as 



" de\ ' 


1 


de% 


~~ E 



1 



- a 
■v 1 



da 1 
da 3 



(22) 



The additional parameter a is included here to take 
into account the anisotropy. When a = 0, we recover 
Hooke's law of Eq. (21). Eq. (17) is calculated from Eq. 
(22) by performing the transformation in the coordinates 
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of the volumetric strain de = —de\ — dti and deviatoric 
strain d^y = de2 — dei, and the pressure p = (01 +02)12 
and the deviatoric stress q = {o\ — <72)/2. One obtains: 



" de e ' 


2 




~ E 



1 — v —a 
—a 1 + v 





dp 




dq 



(23) 



In the isotropic case a = this matrix is diagonal. 
The inverse of the diagonal terms are the bulk modulus 
K = E/2(l - v) and the shear modulus G = E/2(l + v). 
The anisotropy a^O couples the compression mode with 
the shear deformation such that the compression of the 
sample will produce deviatoric deformation. This cou- 
pling can be observed from the inspection of the elastic 
part of the strain envelope responses de e (0) as shown in 
Fig. 6. For stress values such as q/p < 0.4 the stress 
envelope responses collapse on to the same ellipse. This 
response can be described by Eq. (23) taking a = 0. 
For stress values satisfying q/p > 0.4 there is a coupling 
between compression and shear deformations and it is 
necessary to take a ^ in Eq. (23). 



C. Stiffness and fabric 

Comparing the calculation of the elastic response in 
Fig. 6 with the anisotropy of the contact network shown 
in Fig. 5, a certain correlation is evident between the 
stiffness tensor and the fabric coefficients of Eq. (19). We 
observe that Hooke's law is valid in the regime q/p < 0.4 
where the contact network is isotropic. Moreover, we 
observe that the opening of the contacts, whose branch 
vectors are almost perpendicular to the direction of the 
load, produces a reduction of the stiffness in this direc- 
tion. This results in an anisotropic elasticity. 
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+ 9=180° 
9=270° 




FIG. 6. Elastic strain envelope responses de e (8). They are 
calculated for a pressure p = 160fcPa and taking deviatoric 
stresses with q = O.Op (inner), O.lp, ...,0.7p (outer). 

We are going to find a simple relation between the 
orientational distribution of the contacts and the param- 
eters of the stiffness. These three parameters are calcu- 
lated from the elastic response by the introduction of the 
quadratic form of D^ 1 : 



R(9) 



= a T D- 1 a 



dp de e + dq d'y" 
dp 2 + dq 2 



(24) 



Here a T is the transpose of a, which is defined in Eq. 
(15). This function can be directly obtained from the 
elastic part of the strain response de e (0). On the other 
hand, replacing Eq. (23) in Eq. (24) one can express R 
in terms of the parameters of the stiffness tensor: 



R{6) = — [1 - v cos(26») - a sin(20)] . (25) 

Using this equation, the parameters E, v and a are 
evaluated from the Fourier coefficients of R: 



1 



1 



E 4tt 



R{0)d0, 



!/ = -—/ R(6)co5(26)d6, 
2tt Jo 

a = / R(9) sin(20)d0. 

27r Jo 



(26) 
(27) 
(28) 



Figs. 7, 8 and 9 show the results of the calculation 
of the Young modulus E, the Poisson ratio v and the 
anisotropy factor a. Below the limit of isotropy, Hooke's 
law can be applied: E « Eq, v ks v§ and a ~ 0. On 
the other hand, above the limit of isotropy a reduction 
of the Young modulus is found, along with an increase 
of the Poisson ratio and the anisotropy factor. In order 
to evaluate the dependence of these parameters on the 
fabric coefficients from Eq. (19), we introduce an in- 
ternal variable measuring the degree of anisotropy. This 
variable is denoted by a and is defined as the percentage 
change of the total number of contacts. 



N Q -N 
a = = 1 



(29) 



x 10" 



where ao is defined in Eq. (19). The dependence of the 
parameters of the stiffness tensor on the internal variable 
a is evaluated by developing the Taylor series around 
a = 0: 

E{a) = E(0) + E'(0)a + O (a 2 ) , 

a(a) = a(0) + a' (0)a + O (a 2 ) , (30) 

v{a) = u(0) + z/(0)a + v"{0)a 2 + O (a 3 ) . 
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FIG. 7. Young modulus. The solid line is the linear ap- 
proximation of E(d). See Eq. (31). 



FIG. 9. Anisotropy parameter. The dashed line is the lin- 
ear approximation of a(d). See Eq. (31). 
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FIG. 8. Poisson's ratio. The dashed line is the quadratic 
approximation of v(d). See Eq. (31). 



The coefficients of this expansion are caicuiated from 
the best fit of those expansions. Figs. 7 and 9 show 
that the linear approximation is good enough to repro- 
duce the Young modulus and the anisotropy factor. The 
fit of the Poisson ratio, is shown in Fig. 8. The fitting 
with only one internal variable requires the inclusion of a 
quadratic approximation. To obtain more accurate rela- 
tions, it may be necessary to introduce a more complex 
dependence on the fabric coefficients of Eq. (19). 



V. PLASTIC DEFORMATIONS 



In the elasto-plastic models of soils the plastic defor- 
mation is calculated by introducing a certain number of 
hypothetical surfaces [25-27,18]. In the Drucker-Prager 
models, the so-called plastic flow rule is calculated from 
the yield surface and the plastic potential [25,26,18]. In 
the bounding surface plasticity, it is calculated from the 
loading surface and bounding surfaces [27] We will see 
that it is possible to calculate the parameters of the flow 
rule of plasticity, the flow direction, the yield direction 
and the modulus of plasticity, directly from the stress en- 
velope response de p (9) without introducing such abstract 
surfaces. 
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A. Flow rule 

In Fig. 3 we found that the plastic envelope response 
lies almost on a straight line, as is predicted by the 
Drucker-Prager theory [28]. This motivates us to de- 
fine the parameters describing the plasticity in the same 
way as this theory: i.e. the yield direction <f>, the flow 
direction ip, and the plastic modulus h. 

The yield direction is given by the incremental stress 
direction <j> with maximal plastic deformation 



\d<?{<j>)\ = max \d?(6)\, 



(31) 



The flow direction is defined from the orientation of the 
plastic response at its maximum value 



ip = atan2(d<y p , de p ) \ 9=<t> (32) 

Here atan2(y, x) is the four quadrant inverse tangent 
of the real parts of the elements of x and y.( — it <= 
atan2(y, x) <— it). The plastic modulus is obtained from 
the modulus of the maximal plastic response. 

h \da\ 



(33) 



The incremental plastic response can be expressed in 
terms of these quantities. Let us define the unitary vec- 
tors ip and ip . The first one is oriented in the direction 
of ip and the second one is the rotation of %j) of 90°. The 
plastic strain is written as: 
I 



de"(0) 



(34) 



where the plastic components K\{6) and K>2{&) are given 
by 



Kl (6) = hide" ■ ip) 
K 2 (0) = h(de p ■■^ ± ). 



(35) 



These functions are calculated from the resulting plas- 
tic response taking different stress values. The results are 
shown in Fig. 10. We found that the functions Ki(6 — ip) 
collapse on to the same curve for all the stress states. 
This curve fits well to a cosine function, truncated to 
zero for the negative values. The profile k 2 depends on 
the stress ratio we take. This dependency is difficult to 
evaluate, because the values of this function are of the 
same order as the statistical fluctuations. In order to 
simplify the description of the plastic response, the fol- 
lowing approximation is made: 



K 2 (0) « (COS(0 -$)) 



(36) 



where (x) = xQ(x), being Q(x) the Heaviside step 
function. Now, the flow rule results from the substitu- 
tion of Eqs.(34) and (36) into Eq. (18): 



d£ p (6) = J(6)da 



da) 



(37) 



Although we have neither introduced yield functions 
nor plastic potentials, we recover the same structure 
of the plastic deformation obtained from the Drucker- 
Prager analysis [28]. This result suggests the possibil- 
ity to measure such surfaces directly from the envelope 
responses without need of an a-priori hypothesis about 
these surfaces. The next step is to verify the validity 
of the Drucker-Prager normality postulate, which states 
that the yield function must coincide with the plastic 
potential function [28]. 



B. Normality postulate 




-150 ^100 ^50 50 100 150 
9-<|> 



FIG. 10. Plastic components ki(9) (circles) and K2(9) 
(dots). The results for different stress values have been super- 
posed. The solid line represents the truncated cosine function. 



The Drucker normality postulate was established to 
describe the plasticity in metals [28]. The question nat- 
urally arises as to its validity for the plastic deformation 
for soils. With this aim, the yield direction and the flow 
direction have been calculated for different stress states. 
The results prove that both angles are quite different, as 
shown in Fig. 11. A large amount of experimental evi- 
dence has also indicated a clear deviation from Drucker 's 
normality postulate [29]. 

It is not surprising that the Drucker postulate, which 
has been established for metal plasticity, is not fulfilled in 
the deformation of granular materials. Indeed, the prin- 
cipal mechanism of plasticity in granular materials is the 
rearrangement of the grains by the sliding contacts. This 
is not the case of micro-structural changes in the met- 
als, where there is no frictional resistance [30]. On the 
other hand, the sliding between the grains can be well 
handled in the discrete element formulation, which more 
adequately describes the soil deformation. 
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FIG. 11. The flow direction and the yield direction of the 
plastic response. Solid lines represent a linear fit. 

The fact that the Drucker postulate is not fulfilled in 
the deformation of the granular materials has led to the 
so-called non-associated theory of plasticity [18]. This 
theory introduces a yield surface defining the yield direc- 
tions and a plastic potential function, which defines the 
direction of the plastic strain. Both, yield surfaces and 
plastic potential function can be calculated from the yield 
and flow direction, which in turn are calculated from the 
strain envelope response using Eqs. (31) and (32). Ac- 
cording to Fig. 11, they follow approximately a linear 
dependence with the stress ratio q/p: 



= (/><> 



if 1 



p 



(38) 



3.3° ± 



The four parameters ipo = 46° ± 0.75°, ip' 
0.6°, 0o = 78.9°±0.2° and (f>' a = 59.1°±0.4° are obtained 
from a linear fit of the data. This linear dependence with 
the stress ratio has been shown to fit well with the ex- 
perimental data in triaxial [31] and biaxial [32] tests on 
sand. In fact, this implies that the plastic potential func- 
tion and the yield surfaces have the same shape, inde- 
pendent on the stress level. This is a basic assumption 
from several constitutive models [25,26]. 

From Eq. (38), one can see that there is a transition 
from contractancy to dilatancy around q/p = 0.5. This 
transition is typically observed in dense sand under biax- 
ial loading [26] . A consequence of this linear dependency 
is that ip ^ when q = 0. This implies the existence 
of deviatoric plastic strain when the sample is initially 
under isotropic loading conditions, which has also been 
predicted in the original Cam-clay model [25]. 



The existence of deviatoric plastic deformation under 
extremely small loading appears to be in contradiction to 
the fact that the contact network remains isotropic below 
a certain stress ratio. This matter has also been discussed 
by Nova [26], who introduced some modifications in the 
Cam-clay model in order to satisfy the isotropic condi- 
tion [26]. However, we are going to show from a micro- 
mechanical inspection that the oricntational distribution 
of the sliding contacts departs from the isotropy for ex- 
tremely small deviatoric loadings. 



C. Plasticity & sliding contacts 



Under small deformations, the individual grains of a re- 
alistic soil behave approximately rigidly, and the plastic 
deformation of the assembly is due principally to sliding 
contacts (eventually there is fragmentation of the grains, 
which is not going to be taken into account here). A 
complete understanding of soil plasticity is possible, in 
principle, by the investigation of the micro-mechanical 
arrangement between the grains. We present here some 
observations about the anisotropy induced by loading in 
the subnetwork of the sliding contacts. This investiga- 
tion will be useful to understand some features of plastic 
deformation. 

The sliding condition at the contacts is given by \ft\ — 
/ifn, where f n and ft are the normal and tangential com- 
ponents of the contact force, and /x is the friction co- 
efficient. When the sample is isotropically compressed, 
we observe a significant number of contacts reaching the 
sliding conditions. If the sample has not been previously 
sheared, the distribution of the orientation of the branch 
vectors of all the sliding contacts is isotropic. 

This isotropy, however, is broken when the sample is 
subjected to the slightest deviatoric strain. In effect, at 
the very beginning of the loading, most of the sliding 
contacts whose branch vector is oriented nearly parallel 
to the direction of the loading leave the sliding condition. 
The anisotropy of the sliding contacts is investigated by 
introducing the polar function fl B (ip), where £l s ((p)Aip is 
the number of sliding contacts per particle whose branch 
vector is oriented between if and ip + Aip. This can be 
approximated by a truncated Fourier expansion: 



^ S ( < P) ~ 7T^ [ c + Cl COs(2lf) + C2 cos(4<p)] . 



2vr 



(39) 
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FIG. 12. Left: force distribution for the stress ratios 
q/p = 0.1,0.35,0.65. Here ft and f n are the tangential and 
normal components of the force. They are normalized by the 
mean value of f n . Right: orientational distribution of the con- 
tacts fi(vs) (outer) and of the sliding contacts fi s (</?) (inner). 
tp represents the orientation of the branch vector. 



Fig. 12 shows the orientational distribution of sliding 
contacts for different stress ratios. For low stress ratios, 
the branch vectors I of the sliding contacts are oriented 
nearly perpendicular to the loading direction. Sliding oc- 
curs perpendicular to £, so in this case it must be nearly 
parallel to the loading direction. Then, the plastic defor- 
mation must be such as de?, <C de\, so Eq. (32) yields a 
flow direction of %/) rs 45°, in agreement with Eq. (38). 

Increasing the deviatoric strain results in an increase 
of the number of the sliding contacts. The average of 
the orientations of the branch vectors with respect to the 
load direction decreases with the stress ratio, which in 
turn results in a change of the orientation of the plas- 
tic flow. Close to the failure, some of the sliding contacts 
whose branch vectors are nearly parallel to the loading di- 
rection open, giving rise to a butterfly shape distribution, 
as shown in Fig. 12. In this case, the mean value of the 
orientation of the branch vector with respect to the direc- 
tion of the loading is around tp = 38°, which means that 
the sliding between the grains occurs principally around 



52° with respect to the vertical. This provides a crude 
estimate of the ratio between the principal components 
of the plastic deformation as de?, ~ — de\ tan(52°). Ac- 
cording to Eq. (32) this gives an angle of dilatancy of 
'tjj = atan2(d / y p , de p ) ks 97°. This crude approximation 
is reasonably close to the angle of dilatancy of 103.4° 
calculated from Eq. (38). 

A fairly close correlation between the orientation of 
the sliding contacts and the angle of dilatancy has also 
been reported by Calvetti et al. [21] using molecular dy- 
namic simulations in triaxial tests. This correlation sug- 
gests that the plastic deformation of soils can be micro- 
mechanically described by the introduction of the fab- 
ric constants Ci of the equation (39) in the constitutive 
equations. This investigation would lead to some exten- 
sions of the fabric tensor capturing the non-associativity 
of plastic deformation. 
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FIG. 13. Fabric coefficients of the distribution of the con- 
tact normal vectors. They are defined in Eq. (39). 



D. Plastic modulus 

The plastic modulus h defined in Eq. (33) is related 
to the incremental plastic strain in the same way as 
the Young modulus is related to the incremental elas- 
tic strain. Thus, just as we related the Young modulus 
to the average coordination number of the polygons, it 
is reasonable to connect h to the fraction of sliding con- 
tacts n s = N s /N. Here N and N s are the total number 
of contacts and the number of sliding contacts. 

Fig. 14 shows the relation between the hardening and 
the fraction of the sliding contacts. The results can be 
fitted to an exponential relation 

h = h exp(-n s /n ) (40) 
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FIG. 14. Inverse hardening modulus h (MPa) versus frac- 
tion of sliding contacts n s . The values are taken from 
q = 0,0. lp, ...0.7p with different pressures. The lowest value 
of n s corresponds to q — 0. 

Where h = 80.0 ± QAGPa and n = 0.066 ± 0.002. 
This exponential dependence contrasts with the linear 
relation between the Young modulus and the number of 
contacts obtained in Sec. IV. From this comparison, it 
follows that when the number of contacts is such that 
n s > no, the deformation is not homogeneous, but is 
principally concentrated more and more around the slid- 
ing contacts as their number increases. 

The above results suggest that it is possible to estab- 
lish a dependency of the flow rule on the anisotropy of the 
subnetwork of the sliding contacts. This relation is more 
appropriate than just an explicit relation between the 
flow rule and the stress, which probes to be only valid 
in the case of monotonic loading [31]. Since the stress 
can be expressed in terms of micro-mechanical variables, 
branch vectors and contact forces, the identification of 
those internal variables measuring anisotropy and force 
distribution would provide a more general description of 
the dependence of the flow rule on the history of the de- 
formation. 

VI. CONCLUDING REMARKS 

The effect of the anisotropy of the contact network on 
the elasto-plastic response of a Voronoi tessellated sam- 
ple of polygons has been investigated. The most salient 
aspects of this anisotropy are summarized as follows: 

• The incremental elastic response has a centered el- 
lipse as an envelope response. Below the stress ra- 
tio q/p < 0.4, this response can be described by 
the two material parameters of Hooke's law of elas- 
ticity: the Young modulus and the Poisson ratio. 



Above this stress ratio there is a dependence of the 
stiffness on the stress ratio, which can be connected 
to the anisotropy induced in the contact network 
during loading. We should state that this result 
might be dependent on the preparation procedure. 
In particular, samples with void ratio different from 
zero show a smooth transition to the anisotropy, 
which requires further studies. 

• The plastic envelope responses lie almost on the 
straight line defining the plastic flow direction if>. 
The yield direction tjj and the plastic modulus h 
have also been calculated directly from the plastic 
response. The flow direction and yield direction 
depend on the stress ratio. In particular, the plas- 
tic flow for zero stress ratio has a non-zero devia- 
toric component suggesting an anisotropy induced 
for extremely small deviatoric strains. We found 
that this effect comes from the fact that the sliding 
contacts depart from anisotropy when the sample is 
subjected to the smallest deviatoric deformations. 
We have also shown a correlation between the mean 
orientation of the sliding contacts and the flow di- 
rection of the plastic deformations. 

• In the investigation of the connection between the 
plastic deformation and the number of sliding con- 
tacts, we found that the plastic modulus h decays 
exponentially as the fraction of sliding contacts in- 
creases. This contrasts with the linear decrease 
of the Young modulus E with the increase of the 
number of open contacts, suggesting that the de- 
formation of the granular assembly is concentrated 
around the sliding contacts. 

Since the mechanical response of the granular sample 
is represented by a collective response of all the con- 
tacts, it is expected that the constitutive relation can 
be completely characterized by the inclusion of some in- 
ternal variables, containing the information about the 
micro-structural arrangements between the grains. We 
have introduced some internal variables measuring the 
anisotropy of the contact force network. The fabric coef- 
ficients cii, measuring the anisotropy of the network of all 
the contacts, prove to be connected with the anisotropic 
stiffness. On the other hand, the fabric coefficients Cj, 
measuring the anisotropy of the sliding contacts, are re- 
lated to the plasticity features of the granular materials. 

Future work should be oriented towards the elabora- 
tion of a theoretical framework connecting the consti- 
tutive relation to these fabric coefficients. To provide 
a complete micro-mechanically based description of the 
elasto-plastic features, the evolution equations of these 
internal variables must be included in this formalism. 
This theory would be an extension of the ideas which 
have been proposed to relate the fabric tensor to the con- 
stitutive relation of granular materials [33,2,10]. 
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